The contact process in heterogeneous and weakly-disordered systems 
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The critical behavior of the contact process (CP) in heterogeneous periodic and weakly-disordered 
environments is investigated using the supercritical series expansion and Monte Carlo (MC) simu- 
lations. Phase-separation lines and critical exponents j3 (from series expansion) and n (from MC 
simulations) are calculated. A general analytical expression for the locus of critical points is sug- 
gested for the weak-disorder limit and confirmed by the series expansion analysis and the MC 
simulations. Our results for the critical exponents show that the CP in heterogeneous environments 
remains in the directed percolation (DP) universality class, while for environments with quenched 
disorder, the data are compatible with the scenario of continuously changing critical exponents. 

PACS numbers: 05.70.Ln,64.60.Ht,02.50.Ey,87.18.Bb 
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Phase transitions in non-equilibrium statistical me- 
chanics have long been a field of interest, and univer- 
sality classes similar to those in equilibrium have been 
identified, with the DP class being one of the most promi- 
nent ones. The CP EI> a susceptible-infected-susceptible 
model for the spread of epidemics, belongs to this uni- 
versality class and has become one of its archetypical 
models. Recent years have seen much activity sparked 
by the question whether the DP class is robust with re- 
spect to the introduction of quenched spatial disorder or 
not. This has been seen as an important question both 
from a fundamental and experimental viewpoint, as it 
has been suggested that the lack of experimental obser- 
vations of the DP critical behaviour might be a result of 
the presence of disorder in real- world systems |2j. 

One of the foremost arguments that disorder changes 
the critical behaviour of the CP is that it violates the 
Harris criterion 0, 0] for all dimensions d < 4. This cri- 
terion states that a critical point is stable with respect 
to disorder if dv±_ > 2 where v± is the critical exponent 
associated with the spatial correlation length. So far, 
all studies carried out on the disordered CP have pro- 
vided supporting evidence for a change in univers ality 
with the introduction of disorder H M, II H H EH Elf- 
However, it is not entirely clear how the critical ex- 
ponents change with disorder. In the strong-disorder 
limit, Hooybhergs et al. 0, Ell have demonstrated that 
the CP changes to the universality class of the random 
transverse-field Ising model with activated scaling char- 
acterized by known scaling exponents. Recent MC simu- 
lations EH suggest that the activated scaling holds for an 
arbitrary degree of disorder meaning that an introduction 
of even weak disorder forces the abrupt change of critical 
exponents from known values for the homogeneous CP 
to those of the infinite-randomness fixed point (IRFP). 
This contradicts the findings of other authors who 
showed, both using MC simulations and density-matrix 



renormalization-group (RG) analysis, that there is an in- 
termediate disorder regime with continuously varying ex- 
ponents. Unconventional critical behaviour produced by 
quenched randomness is supported by a field-theoretical 
analysis in which only runaway solutions in the RG 
equations were found. 

The subject of this paper is to investigate the Id CP 
in heterogeneous periodic systems (e.g. a regular binary 
chain) and in systems with weak disorder (e.g. a binary 
chain with randomly placed species characterized by pa- 
rameters close in value) . In order to achieve this goal, we 
employ the supercritical series expansions [T3 . E5j | and 
MC simulations. We also suggest a simple analytical ex- 
pression for the locus of critical points in the rate-space 
phase diagram which is in very good agreement with se- 
ries expansion and MC simulation results both for het- 
erogeneous periodic and weakly-disordered systems. Our 
main findings demonstrate that the CP in heterogeneous 
periodic systems belongs to the DP universality class 
with the scaling exponents coinciding with those for ho- 
mogeneous case. For weakly-disordered systems, we can 
state that the introduction of disorder does not force the 
exponents to change to the values of the IRFP but rather 
causes their continuous change with disorder strength. 

The CP is usually defined on a hypercubic lattice of 
nodes which can be either empty (susceptible) or occu- 
pied (infected). The infection occurs via contacts be- 
tween Z nearest nodes i and j with the rate Xij/Z. An 
infected node i can recover to susceptible one with the 
rate \Xi. The time scale is defined by setting all \j = 1 
(for simplicity, there is no disorder in transmission rates) 
and, for concreteness, we consider only binary systems 
with two types of nodes, A and B, characterized by the 
recovery rates, /ia and /ib, respectively, which are dis- 
tributed according to bimodal distribution in the disor- 
dered system, p(fii) = (1 - q)8(jj,i - ha) + qS(p-i - hb), 
with q being the concentration of nodes B. 
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In the homogeneous case (q = 0), the CP undergoes a 
non-equilibrium phase transition between active and ab- 
sorbing states [14| if the recovery rates become greater 
than a critical value, /i > /i c ~ 0.303228 At crit- 

icality, the number of infected sites, N\„f(t), scales with 
time as, N mf (t) oc t r i, with rj = 0.313686 15]. Close to 
criticality, the survival probability, P^, and the time and 
space correlation lengths, £| and £x> respectively, also ex- 
hibit typical critical behaviour, P^ oc A*, £|| oc |A| _!y n 
and £j_ oc |A| _l/i where A = /z c — \x and the universal 
exponents are (3 ~ 0.2769 v\\ = 1.733825(25) and 
u x = 1.096844(14) 

In the heterogeneous and disordered systems, a similar 
transition occurs 0. Below, we address two questions: 

(i) how heterogeneity and disorder influence the universal 
properties (namely, scaling exponents) at criticality and 

(ii) what the locus of critical points is in the rate space 
of the CP in such systems. We start with simple argu- 
ments about a possible way of evaluating an analytical 
expression for the critical line separating the active and 
absorbing states. 

Let us consider a system of N (N — > oo) nodes char- 
acterized by random recovery rates, jli = /ii//i c . The 
locus of critical points in the space of recovery rates can 
be written in a general form as a solution of the follow- 
ing equation, F(M{) = (for i = 1,2, . . . , N), where 
Mi = ln/2j. The function F{Mi) is invariant under the 
exchange of any two arguments and it obeys the prop- 
erty f(0) = for the homogeneous case. Assuming that 
this function is analytic around the homogeneous crit- 
ical point we can expand it in a Taylor series around 
this point, F'{0)Y, l M t + O(Mf) = 0, where we have 
used the symmetry of the function F leading to all first 
derivatives being equal at the stationary point. In fact, 
this expansion is equivalent to the expansion in the mo- 
ments u n — E[M n ], where E[-] denotes the expectation 
value. Leaving only the first order in the above expan- 
sion we end up with the following approximate equation 
for the locus of critical points, 

£[ln/x]~0, (1) 

which is valid around the homogeneous critical point. 
The choice of arguments Mj of function F is motivated 
by a similar choice for the RG analysis of both the CP 
and random transverse-field Ising model 0,^3- Below, 
we demonstrate that Eq. (JTJ indeed describes well the 
phase-separation lines in all studied heterogeneous and 
disordered systems. 

In order to support Eq. (JTJ for the locus of critical 
points and also to investigate the scaling behaviour at 
criticality, we have used the perturbative supercritical 
series expansion 0, for the survival probability, ■ 
Following the formalism developed in Refs. [lSL Il9l l2fl|. 
the state of the system is described by the state vec- 
tor \P(t)) = J2fa\ P{{ a }^) whose time evolution 



is governed by the master equation d t \P(t)) = L\P(t)). 
Here L is the generator of the Markov process that 
contains the transition rates between the different mi- 
crostates of the system, \{cr}). As usual, L = fiW + V , 
is split into a part that destroys particles, fiW, and a 
part that creates particles, V, which in the systems un- 
der consideration take the forms 

fiW = ^(l-atja, (2) 

i 

V = XI ^ C 1 ^ a i ) a J + «J+i a i+i) . ( 3 ) 

i 

where a] and a, are hard-core bosonic creation and an- 
nihilation operators, respectively. The Laplace trans- 
form of \P(t)), \P(s)) = (s - nW - Vy^PiO)), is 
then expanded in /Zj, yielding the supercritical expansion 
for the survival probability P 00 (/i J 4, hb) = hm s ^ (l — 
s(0|P(s)}), where |0) is the absorbing state. For the 
analysis of this multi-variable survival probability (cf. 
Ref. H^l), we employ a numerical scheme similar to the 
Nested Pade Approximation |22t |23| . In order to investi- 
gate the critical behaviour, we consider the meromorphic 
function d fJ , A \nP 00 (iiA,S) (with 5 — \ib — Ha), whose 
first pole on the positive real axis is the critical point 
of the model and the residue at that pole is the critical 
exponent [3. To improve estimates of these poles from 
the finite scries expansion, the following multivariable 
rational-approximant scheme is used: for a given expan- 
sion of Poo in two variables, ha and 5, up to an even 
(odd) order N, the Pade approximants [n, n] ([n, n + 1]) 
in 5 of the coefficients of the terms /i^~ 1-2n (/^ -2n ) in 
the series d )lA lnP 00 (/ZA, 6) are formed, followed by the 
construction of the Pade approximant [N/2 — 1, N/2] 
([(N - l)/2, (N - l)/2]) in fi A . In order to estimate the 
stability of the poles and residues found, several Pade ap- 
proximants in fx a (e.g. the approximants from [N/2 — 1, 
N/2] down to [N/2 - 2, N/2 - 1] for even orders of N) 
were constructed and averaged over. 

Using this scheme and starting from a single seed in 
the series expansions up to order N — 24, we calcu- 
lated the locus of critical points and the critical ex- 
ponents (3 for three heterogeneous lattices, AB, AAB, 
and AABB, and for disordered systems whose recov- 
ery rates are drawn from the bimodal distribution men- 
tioned above (see Figs. Q] and HJ). Fig. Qfa) demonstrates 
that around the homogeneous critical point (fx c , (i c ), the 
phase-separation lines between the active and absorbing 
state are indeed very well described by Eq. (JJJ. This 
is also confirmed by single-seed MC simulations based 
on the random-sequential algorithm [24| (see Table P) 
- the deviations of the MC results from predictions of 
Eq. JJJ and series expansion data for critical line are less 
than 1%. The critical values of [Ib for fixed values of 
fiA were found by analysis of the power-law behavior 
of Ni n [(t) with averaging over 10 6 MC runs. Fig. ^b) 
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FIG. 1: (Color online) Periodic Id lattices, AB (o), AAB 
(o), and AABB (+): (a) critical points obtained by series 
expansions in comparison with analytical prediction for the 

critical line fi c = (uA~ q lA), 9 = 1/2 ( — ) and q = 1/3 ( ) 

for fi c = 0.303228 Qjj]; (b) critical exponent (3 in compari- 
son with series expansion value (3 = 0.2769 ( ) 13] for the 

homogeneous case. 



TABLE I: The critical values of /is obtained in MC simula- 
tions (second column) and calculated according to Eq. Q 
(third column) together with the MC critical exponents 
(fourth column) for fixed values of {ia (first column) in the 
heterogeneous systems AB and AAB. 



HA 



Hb (MC) (pred.) 



'/ 



AB 



AAB 



0.2750 
0.2500 
0.2250 
0.2750 
0.2500 
0.2250 



0.3344 ± 0.0001 
0.3681 ± 0.0001 
0.4094 ± 0.0001 
0.3689 ± 0.0001 
0.4475 ± 0.0001 
0.5546 ± 0.0001 



0.3343 
0.3678 
0.4087 
0.3687 
0.4461 
0.5507 



0.313 ± 0.006 
0.313 ± 0.003 
0.313 ± 0.003 
0.313 ± 0.002 
0.314 ± 0.001 
0.314 ± 0.001 



and Tab. |U also show that the critical exponents j3 and r\ 
practically do not change from the values for the homo- 
geneous CP, thus confirming that the CP in the heteroge- 
neous lattices belongs to the same universality class, DP, 
as the homogeneous one. The systematic deviations of 
the calculated critical rates from the theoretical predic- 
tion increases with the distance from the homogeneous 
point thus reflecting the restricted range of applicabil- 
ity of Eq. (JJJ. Some irregular fluctuations both in hb 
and (3 are probably due to poor convergence of the series 
expansions. 

A similar series expansion analysis has been performed 
for disordered systems with a configurational averaging 
of the survival probability, (Poo). We were able to per- 
form the complete numerical averaging over all 2 2N ~ 1 
configurations for N < 12. The configurational aver- 
aging for series expansions to higher orders, N = 19, 



FIG. 2: (Color online) The phase diagram (a) and scaling 
exponent /3((j,a), (b) and inset of (a), for disordered lattices 
with various degrees of disorder: q = 0.04 (o), q = 0.3 (A) 
and q = 0.5 (□) obtained by series expansion. The lines 
represent the theoretical prediction according to Eq. Q , /i c = 

{^ A ~ q (1%), for q = 0.04 (— ), q = 0.3 ( ) and q = 0.5 (• • • ) 

with )i c — 0.303228 fl3|. The triangle in (a) marks the region 
for which the MC simulations shown in Fig.[3]have been run. 
The dashed lines in (b) and the inset of (a) show the value 
of /3 for the homogeneous case, (5c = 0.2769 (---) 13]. The 
arrow in (b) marks the homogeneous critical point. 



has been done approximately. In the calculation of 

(•Poo) = EtoELo( c "-»>/ i B/ 1 i" ra > at each term of or- 
der M < N, we only included realizations with num- 
ber i < %m of "impurity" -sites B: e.g. we have chosen 
iN—n = n + 2 for series expansions up to order N = 19, 
so that for M = 18 only realizations with up to three 
impurities contributed to {Poo)- Dropping these disor- 
der realizations from the configurational average incurs 
less error the smaller the impurity concentration q is. 
Assuming that the coefficients cum ~ cm (with cm be- 
ing the coefficient of the same order in the homogeneous 
case), all the realizations with i impurities contribute a 
term ( 2 ^^ _1 )(1 — q) 2AI ~ 1 ~' l q' 1 to the configurational av- 
erage, (cMm)- Then, it can easily be shown that for 
1 = 9max = 0.04 all terms with i > — 2 are smaller 
than the terms with i < iff. The validity of this approxi- 
mation has been confirmed by testing it against the exact 
results for N < 12. 

The results for fully and partially averaged survival 
probabilities in disordered systems are shown in Fig. [21 
The phase-separation lines have been obtained for arbi- 
trary impurity concentration for the fully averaged Poo 
expanded up to order N — 12 and two of them for q = 0.5 
(squares and dotted line) and q = 0.3 (triangles and 
dot-dashed line) are displayed in Fig. Eta). High-order 
series expansions (N = 19) have been calculated only 
for low impurity concentrations, q < <7 ma x = 0.04 (see 
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the circles and solid line in Fig. Ufa)). Again, the poles 
of dfj, A In Pqo (//a, S) agree very well with the theoretical 
prediction given by Eq. QJ. 

The residues of the poles (exponents 8) for different 
points on the critical line are shown in Fig. |2[b) and in 
the inset in Fig. Ufa). The value of 8 reaches a mini- 
mum, ft m in, located approximately around the homoge- 
neous critical point (fx c , fJ, c ) with 8 m in being rather close 
to the value of the homogeneous critical exponent, 8 C , 
with |/3 min - 8 C \//3 C ~ 0.5% for N = 12 and 0.21% for 
N = 19 (see the inset in Fig. Ufa)) thus confirming that 
the value of the exponent is much more sensitive to N 
than the critical rates. Away from the critical point, the 
value of 8, first, monotonically increases and then starts 
to fluctuate due to a high sensitivity to the value of \xb, 
the estimates of which loose precision due to poor con- 
vergence of the series in this range. The results for 8 are 
in reasonable agreement with findings in Refs. 0, Q, 0] 
where continuously varying critical exponents were seen 
in MC simulations and density-matrix renormalization- 
group studies of the random CP. Unfortunately, the er- 
rors in the exponents that are shown in Fig. 0b) and the 
inset of Fig. 0a) are at least of the order of \8(fi c ) — 0c\ 
and the monotonic growth of the exponents found in the 
series expansions can be questioned. However, our results 
are certainly not consistent with the scenario presented in 
[lH according to which the weakly-disordered CP belongs 
to the same universality class as random transverse-field 
Ising model with 8 = 0.38197. 

The results for the disorder case have been supported 
by MC simulations (see Fig. 0). Due to the longrelax- 
ation times of the disordered CP (cf. Refs. 0,0, El) 
we have focused only on one point in the rate space with 
q = 0.5 and [Ia = 0.25. The results of the simulations up 
to 10 7 time steps are shown in Fig.0for three values of hb 
around criticality in double-log scales of N[ n f vs t as well 
as iVi n f vs In t (see the inset in Fig. 0) to allow for both 
conventional and activated scaling [llj. The MC criti- 
cal value of hb — 0.368 (with the error being less than 
0.005) obtained by conventional double-log scaling anal- 
ysis is certainly very close to the series expansion value 
(fj, B ~ 0.369 for N = 12, see the triangle in Fig. 01. The 
value of the dynamical exponent at this point is found to 
be r\ ~ 0.388 which is in favour of the scenario suggesting 
scaling exponents varying continuously with disorder. 

In conclusion, we have investigated the CP in hetero- 
geneous and disordered Id systems in the limit of weak 
disorder by means of the series expansions and MC simu- 
lations. We have demonstrated that the CP in heteroge- 
neous Id lattices stays in the DP universality class. For 
disordered environment, our results suggest that disor- 
der continuously changes the scaling exponents. A sim- 
ple analytical formula for the phase-separation line has 
been suggested and proved (numerically) to be valid in 
the weak-disorder limit. Preliminary investigations of the 
CP in 2d heterogeneous lattices also support the analyt- 
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FIG. 3: (Color online) Number of infected sites vs time ob- 
tained by MC simulations of the disordered CP with q = 0.5 
and /xa = 0.25, averaged over 1500 runs and at least 30 disor- 
der realizations. The main figure is in double-log scale accord- 
ing to the conventional scaling while the inset demonstrates 
the same data in the activated scaling picture [i~H| . The curves 
from top to bottom correspond to 0.3628 (A), 0.3680 (□) and 
0.3728 (O). 



ical predictions for the phase-separation line. 

Acknowledgements - The authors would like to thank 
I. Jensen, B. Simons, J. Stilck, and R. Dickman for 
their helpful correspondence and remarks. The compu- 
tations were mainly performed on the Cambridge Uni- 
versity Condor Grid and the Cambridge-Cranfield High- 
Performance Computer facility. CJN and SVF would like 
to thank the EPSRC and the Cambridge European Trust 
for finacial support. 



* Electronic address: cjn24@cam.ac.uk 
[1] T. Harris, Annals of Probability 2, 969 (1974). 
[2] H. Hinrichsen, Advances in Physics 49, 815 (2000). 
[3] A. B. Harris, J. Phys. C 7, 1671 (1974). 
[4] J. Chayes et al, Phys. Rev. Lett. 57, 2999 (1986). 
[5] A. Noest, Phys. Rev. Lett. 57, 90 (1986). 
[6] A. Moreira and R. Dickman, Phys. Rev. E 54, 
R3090(1996). 

[7] R. Dickman and A. Moreira, Phys. Rev. E 57, 1263 
(1998). 

[8] H. K. Janssen, Phys. Rev. E 55, 6253 (1997). 

[9] J. Hooyberghs et al, Phys. Rev. Lett. 90, 100601 (2003). 
[10] J. Hooyberghs et al, Phys. Rev. E 69, 66140 (2004). 
[11] T. Vojta and M. Dickison, Phys. Rev. E 72, 36126 (2005). 
[12] R. Dickman and I. Jensen, Phys. Rev. Lett. 67, 2391 
(1991). 

[13] I. Jensen and R. Dickman, J. Stat. Phys. 71, 89 (1993). 
[14] T. Liggett, Interacting Particle Systems (1985), 1st ed. 
[15] I. Jensen, J. Phys. A 29, 7013 (1996). 



■5 



[16] M. Bramson et al., Ann. Probab. 19, 960 (1991). 
[17] D. S. Fisher, Phys. Rev. B 51, 6411 (1995). 
[18] M. Doi, J. Phys. A 9, 1479 (1976). 

[19] P. Grassberger and A. De La Torre, Ann. Phys. (N.Y.) 

122, 373 (1979). 
[20] L. Peliti, J. Physique 46, 1469 (1985). 
[21] W. G. Dantas and J. F. Stilck, J. Phys. A: Math. Gen. 



38, 5841 (2005). 
[22] P. Guillaume, J. Comp. Appl. Math. 82, 149 (1997). 
[23] P. Guillaume and A. Huard, J. Comp. Appl. Math. 121, 

197 (2000). 

[24] J. Marro and R. Dickman, Nonequilibmum Phase Tran- 
sitions in Lattice Models (CUP, 1999, 1st ed.). 



